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ABSTRACT 

Context. The electron-cyclotron maser instability is responsible for generation of the auroral kilometric radiation of the Earth and 
similar phenomena at other magnetized planets of the Solar system. Recently discovered radio emission from ultracool dwarfs has 
many similarities with the planetary auroral radio emissions. The in situ measurements in the terrestrial magnetosphere indicate that 
the radiation results from nonthermal electrons with a horseshoe-like distribution. Kinetic simulations of the electron-cyclotron maser 
instability for such distribitions have not been done yet. 

Aims. In this work, we investigate amplification of plasma waves by the horseshoe-like electron distribution as well as relaxation of 
this distribution due to the electron-cyclotron maser instability. We aim to determine the parameters of the generated plasma waves, 
timescales of the relaxation process, and the conversion efficiency of the particle energy into waves. 

Methods. We have developed a kinetic relativistic quasi-linear 2D code for simulating the coevolution of an electron distribution 
and the high-frequency plasma waves. The code includes the processes of wave growth and particle diffusion which are assumed to 
be much faster than other processes (particle injection, etc.). A number of simulations have been performed for different parameter 
sets which seem to be typical for the magnetospheres of ultracool dwarfs (in particular, the plasma frequency is much less than the 
cyclotron one). 

Results. The calculations have shown that the fundamental extraordinary mode dominates strongly. The generated waves have the 
frequency slightly below the electron cyclotron frequency and propagate across the magnetic field. The final intensities of other modes 
are negligible. The conversion efficiency of the electron energy into the extraordinary waves is typically around 10%. Complete 
relaxation of the unstable electron distribution takes much less than a second. 

Conclusions. Energy efficiency of the electron-cyclotron maser instability is more than sufficient to provide the observed intensity of 
radio emission from ultracool dwarfs. On the other hand, the observed light curves of the emission are not related to the properties of 
this instability and reflect, most likely, dynamics of the electron acceleration process and/or geometry of the radiation source. 

Key words, radiation mechanisms: non-thermal - planets and satellites: aurorae - brown dwarfs - radio continuum: stars 
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NJ . et al. I20081 . The emission includes a slow-varying weakly- 

■ The electron-cyclotron maser instability (ECMI) can occur in polarize d component as well as short intense periodic bursts with 

a magnetized plasma with a non-equilibrium electron distnbu- almost m% po^atio^ whose period seems to coincide with 

O tion function (a positive slope in perpendicular velocity is re- the rotational p er i d of the star. While the slow-varying compo- 

. quired) and results in an effective amplification of electromag- nent may be exp lained by the incoherent gyrosynchrotron radia- 

L; ; netic waves at the low harmonics of the elec tron cy clotron fre- tion> the per iodic bursts require a coherent radiation mechanism 

. r-j . quencyJWu & Lee [1979J Melrose & Dulk [1981 Winglee & such as ECMI Measurements of magnetic fields for the stars of 

^ ■ Dulk [1981. In a relatively rarefied plasma (when the electron late spectral class (M4) using phase-resolved spectropolarimetry 

H ; plasma frequency to cyclotron frequency ratio w p /w B « 1), the (Donati et aL 2 006 Morin et al. 120081 indicate that the magnetic 

, dominant mode of astrophysical electron-cyclotron masers is the field has a di p ole -iike structure with the dipole axis close to the 

extraordinary wave with the frequency close to the fundamental rotational axis . Thus the ra dio emission of UCDs seems to be 

cyclotron frequency. The emission has a narrow spectral band- similar to the auroral em i ss ions of the Solar system planets, but 

width, a narrow directivity pattern, and high (nearly 100%) de- with much a stronger magnetic field (not less than 3000 G, to 

gree of circular polarization. ECMI is responsible for generation provide the highest obserV ed emission frequency). Existence of 

of the auroral kilometric radiation (AKR) in the magnetosphere such magnet ic fields at the cool stars (<M9) is also confirmed 

of the Earth and, most likely, the similar phenomena at other via in f rarecl measurements (e.g., Reiners & BasriHED. 
magnetized planets of the Solar system (Zarka 1998, Treumann 

120061 1. ECMI has also been applied to interpret certain types of ECMI has been investigated analytically and using numeri- 

solar and stellar sporadic radio bursts (Melrose & Dulk 119821 C al simula tions i n a number of papers (see, e.g., the review of 

Fleishman & Melnikov[l998}. Treumann 120061 and references therein). The required unstable 

Recently, a number of late M stars and brown dwarfs (termed electron distributions (with a deficiency of particles with low 

ultracool dwarfs, UCDs) have been found to be the sources of transversal velocity) are formed naturally due to particle reflec- 

unexpectedly intense radio emission at the frequencies of a few tion from a magnetic field gradient. Using a linear approxima- 
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tion, it is easy to show that even relatively weak fluxes of accel- 
erated electrons can result in large growth rates, so the waves 
can be amplified from the level of thermal fluctuations up to 
the observed intensities at the distances not exceeding a few 
kilometers (see, e.g., Melrose & Dulk 1982, Bingham & Cairns 
2000 Bingham, Cairns, & Kellett 2001 as well as estimations in 
this article). However, the resulting wave energy should become 
comparable with the particle energy; this results in a consider- 
able particle diffusion on the waves, so the linear approximation 
becomes inapplicable. To interpret the cosmic radio emissions 
by ECMI, we have to use a nonlinear model that takes into ac- 
count relaxation of unstable electron distribution due to interac- 
tion with the excited waves. In particular, only nonlinear models 
can provide us with such an important parameter as the transfor- 
mation coefficient of the particle energy into waves. 

Initially, ECMI was associated with the electron distribu- 
tion of the loss-cone type. Kinetic simulation of relaxation of 
the loss-cone is made in works of Aschwanden (1990) and 
Fleishman & Arzner ( 120001 1. However, in situ measurements in 
the AKR sources have shown that the radiation should be pro- 
duced (at least, in this particular case) by ring-like or horseshoe- 
like distribution (e.g., Delory et al. 1998; Ergun et al. [2000). 
These distributions are formed when electron beams (acceler- 
ated by an electric field) move into stronger magnetic field re- 
gions. Kinetic simulation of relaxation of such distributions have 
not been done yet (although there are some particle-in-cell simu- 
lations; see, e.g., the references in the review of Treumann 2006). 

When simulating relaxation of unstable electron distribu- 
tions, the greatest challenge is spatial movement of the plasma 
waves and particles. Wave propagation out of the region occu- 
pied by the electrons with an unstable distribution naturally lim- 
its the wave amplification. On the other hand, the waves am- 
plified in one part of the radiation source will cause relaxation 
of the electron distribution in other parts. Electron movement in 
an inhomogeneous magnetic field is necessary to form an un- 
stable distribution (of the loss-cone or horseshoe type). To take 
into account all the mentioned factors, one has to create a 3D 
model of the emission source and adjacent regions (e.g., a mag- 
netic loop in the solar or stellar corona) which, in turn, requires 
enormous computational resources. The most popular solution 
of this problem is to neglect the spatial movement of waves and 
particles completely, so that the model is reduced to an initial 
value problem where an arbitrary unstable electron distribution 
is used as the initial condition. We can neglect spatial move- 
ments if the time of wave/particle escaping from the wave am- 
plification region r esc far exceeds the typical diffusion time Tdiff . 
Such approach (a diffusive limit) is used in the above mentioned 
papers on kinetic simulation as well as in most particle-in-cell 
simulations. Obviously, this model is far from reality since it 
cannot explain long-term generation of emission and, in addi- 
tion, requires almost instant formation of an unstable electron 
distribution. Nevertheless, a model considering only the wave 
growth/damping and particle diffusion allows us to (i) investigate 
the qualitative behaviour of relaxation of unstable electron dis- 
tributions, (ii) determine the parameters of the produced waves, 
(iii) estimate the diffusion and relaxation timescales (which is 
necessary to make a conclusion about the validity of the diffu- 
sive limit), and (iv) determine the transformation coefficient of 
the accelerated particles energy into waves. 

Measurements in the AKR sources have shown that the cold 
electron component is almost absent there. Under such con- 
ditions, the wave dispersion is determined mainly by the en- 
ergetic electrons and the relativistic effects become important. 
Relativistic corrections to the dispersion relation result in a de- 



crease in the cutoff frequency of the fast extraordinary mode; in 
addition, the fast and slow branches of the extraordinary mode 
can reconnect to form a single branch (Winglee 119831 [1985; 
Strange way [19851 [1986; Robinson [19861 [19571 Le Queau & 
Louarn [19891 Louarn & Le Queau [19961 1. These effects allow 
the waves generated at the frequencies below the nonrelativis- 
tic cyclotron frequency to escape freely from the source into 
vacuum. Solving the exact relativistic dispersion equation is a 
complicated task. However, the existing studies (e.g., Robinson 
1986, 119871 ) have shown that in a sufficiently hot low-density 
plasma, the wave dispersion becomes like that in vacuum (with 
the refraction index N — > 1 and group speed v gr — > c). Such an 
approximation seems to be valid in the sources of the planetary 
auroral radio emissions. 

In this work, we have developed a kinetic relativistic 2D 
code for simulating the coevolution of electron distributions and 
plasma waves in the diffusive limit. Different wave modes can 
be considered both separately and simultaneously. Electron dis- 
tribution of the horseshoe type is considered as the source of 
plasma oscillations. The calculations are made for the conditions 
that seem to be typical for UCDs, although the results can be eas- 
ily scaled to other emission sources (e.g., the Earth or Jupiter). 
Exact plasma parameters in the magnetospheres of UCDs are un- 
known; in particular, we do not know whether the cold plasma 
component is present or not. Therefore we consider two opposite 
cases: (i) when a low-temperature plasma with the maxwellian 
distribution dominates and (ii) when such a component is absent. 
In the former case, we use the cold plasma dispersion relation; in 
the latter case, the dispersion relation is assumed to be like that 
in vacuum. 

Note that the model used is restricted to the high-frequency 
elecromagnetic/magnetoionic waves and does not consider gen- 
eration of the low-frequency waves (e.g., acoustic ones). In the 
AKR sources, amplitude of the low-frequency waves can reach 
high levels resulting in formation of solitary structures (such as 
electron and ion holes) which, in turn, can affect the generation 
of the radio emission (Pottelette, Treumann, & Berthomier 2001 ; 
Mutel et al. 2006, 2007). However, these effects are beyond the 
scope of this paper. 

The model used is described in Section [2] The initial condi- 
tions of the model (including the electron distribution function) 
are described in Section [3] The simulation results are presented 
in Section [4] and discussed in Section [5] The conclusions are 
drawn in Section [6] The calculation formulae (most of which 
can be found elsewhere) are given in Appendices. 

2. Kinetic equations 

If spatial movement of waves and particles is neglected, the co- 
evolution of the electron distribution and waves in a plasma may 
be described in the general case by the following system of equa- 
tions: 



8W£\k,t) 
dT 



dW (A °(k f) 
d/(p,f) 



= r 



(1) [k,/(p,f)]W^(k,f), 



(i) 
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where w£ is the energy density of oscillations of mode cr in the 
space of wave vectors, cr - 1, . . .,N, y (cr) is the growth rate for 
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/b,Hz 


n b /n 




E b , keV 


a 




4i, 


s 






s 






^"sa t Tmax 


w m /w w 


1 


4x 


10 5 


io- 2 


IO" 2 


10 


60° 


3.12 x IO 2 


2.09 x 


10" 


-1 


1.97 x 


io- 


-1 


65.2 


61.6 


0.118 


2 


4x 


10 5 


1 


IO" 3 


10 


60° 


3.48 x IO 2 


2.53 x 


10" 


-1 


1.94 x 


io- 


-1 


88.2 


67.4 


0.131 


3 


4x 


10 7 


10 2 


io- 2 


10 


60° 


3.12 x 10 4 


1.82 x 


10" 


-3 


1.73 x 


io- 


-3 


56.8 


53.9 


0.119 


4 


4x 


10 7 


1 


io- 3 


10 


60° 


3.48 x IO 2 


1.65 x 


10" 


-3 


1.64 x 


io- 


-3 


57.3 


57.0 


0.132 


5 


4x 


10" 


10~ 4 


io- 2 


10 


60° 


3.12 x 10 4 


1.00 x 


10" 


-3 


1.13 x 


io- 


-3 


31.2 


35.1 


0.116 


6 


4x 


10" 


10~ 4 


io- 1 


10 


60° 


2.79 x IO 5 


1.51 x 


10" 


-4 


2.79 x 


io- 


-4 


42.2 


77.8 


0.107 


7 


4x 


10" 


IO 2 


io- 2 


3 


60° 


5.33 x 10 6 


6.40 x 


10" 


-6 


7.15 x 


io- 


-6 


34.1 


38.1 


0.0841 


8 


4x 


10" 


IO 2 


io- 2 


10 





2.67 x 10 6 


1.80 x 


10" 


-5 


2.09 x 


io- 


-5 


48.0 


55.8 


0.103 


9 


4x 


10" 


IO 2 


io- 2 


10 


60° 


3.12 x 10 6 


1.05 x 


10" 


-5 


1.54 x 


io- 


-5 


32.8 


48.2 


0.120 


10 


4x 


10" 


IO 2 


io- 2 


10 


90° 


2.96 x 10 6 


2.21 x 


10" 


-5 


1.99 x 


io- 


-5 


65.5 


58.8 


0.115 


11 


4x 


10" 


IO 2 


io- 2 


10 


120° 


1.88 x 10 6 


1.91 x 


10" 


-5 


2.46 x 


io- 


-5 


35.8 


46.2 


0.0552 


12 


4x 


10" 


IO" 2 


io- 2 


30 


60° 


1.15 x 10 6 


5.77 x 


10" 


-5 


4.25 x 


io- 


-5 


66.3 


48.8 


0.130 


13 


4x 


10" 




io- 3 


3 


60° 


1.16 x IO 7 


3.12 x 


10" 


-6 


4.39 x 


io- 


-6 


36.2 


51.0 


0.133 


14 


4x 


10" 




io- 3 


10 





2.99 x 10 6 


2.19 x 


10" 


-5 


2.24 x 


io- 


-5 


65.5 


66.9 


0.115 


15 


4x 


10" 




io- 3 


10 


60° 


3.48 x 10 6 


9.45 x 


10" 


-6 


1.44 x 


io- 


-5 


32.3 


50.2 


0.133 


16 


4x 


10" 




io- 3 


10 


90° 


3.32 x 10 6 


9.82 x 


10" 


-6 


1.74 x 


io- 


-5 


32.6 


57.7 


0.129 


17 


4x 


10" 




io- 3 


10 


120° 


2.21 x 10 6 


2.13 x 


10" 


-5 


1.76 x 


io- 


-5 


47.3 


39.0 


0.0621 


18 


4x 


10" 




io- 3 


30 


60° 


1.17 x 10 6 


5.68 x 


10" 


-5 


4.51 x 


io- 


-5 


66.3 


52.7 


0.133 


19 


4x 


10" 




io- 2 


10 


60° 


3.48 x IO 8 


1.65 x 


10" 


-7 


1.64 x 


io- 


-7 


57.3 


56.9 


0.132 



a given mode, k is the wave vector, / is the electron distribution 
function, Dry is the diffusion tensor (describing electron scatter- 
ing on the waves of mode <x), and p is the electron momentum. 
The equations in the system (Q3 are coupled implicitly, since the 
growth rate of plasma waves depends on the electron distribu- 
tion function, while the diffusion tensor depends on the intensity 
and spectral distribution of the plasma waves. The expressions 
for the growth rate and elements of the diffusion tensor are given 
in Appendices IB1 and ICl respectively. 

In this work, we explore coevolution of the electron dis- 
tribution and plasma waves using numerical simulations. The 
distributions /(p) and w£ (k) are defined on regular grids in 
(p, a)- and (oj, ^-spaces, respectively, where a is the electron 
pitch angle, a> is the wave frequency, and 6 is the wave propaga- 
tion direction (with respect to the ambient magnetic field). For 
different modes, different grids are used. In most simulations, 
we have considered two wave modes (e.g., ordinary + extraordi- 
nary). The grid size was chosen to be 60 x 60 data points for all 
considered distributions. Note that in the work of Aschwanden 
( 1990), plasma waves were desribed using an irregular adaptive 
grid where the density of data points in phase space was approxi- 
mately proportional to the initial growth rate. However, our sim- 
ulation have shown that for the ring-like or horseshoe-like elec- 
tron distributions, the region of positive growth in phase space 
can shift noticeably during the process of relaxation. Therefore, 
a regular grid is more suitable, and the considered area in (u>, 6)- 
space has to be wider than the initial region of positive growth. 
The system of equations ([TJ is integrated with respect to time 
using the Gear formulae of fourth order (see AppendixlDlfor de- 
tails of the numerical code). All cyclotron harmonics affecting 
the waves growth/damping and particle diffusion are considered; 
however, as a rule, the effect of the first harmonic is dominant. 

Since the simulated system is closed, energy and particle 
number must be conserved in it (total system energy equals the 
sum of energy of particles and energies of all considered oscil- 
lation modes). Fulfillment of the conservation laws can be con- 



sidered as a test of self-consistency in the model and accuracy 
of the numeric code. In our simulations, at the late stage of the 
relaxation process, the particle number was conserved with the 
relative error < 1 .5 X 10~ 3 , and the total energy of the system was 
conserved with the relative error < 3 x 10~ 3 . These estimations 
correspond to the models without the thermal plasma compo- 
nent; for the models including the thermal plasma, the computa- 
tion accuracy can be even better. 

3. Initial conditions 

In this work, we assume that the initial electron distribution func- 
tion has the form: 



/(p, f)\ t=0 = (n - n b )/o(p) + «b/b(p), 



(2) 



where n is the total electron concentration, «b is the concentra- 
tion of nonthermal electrons, /o is the maxwellian distribution 
function of thermal electrons, and /b is the distribution function 
of nonthermal electrons (both functions /o and /b are assumed to 
be normalized to unity). Both the cases of «b ^ n (thermal com- 
ponent dominates) and «b = n (thermal component is absent) are 
considered. 

As stated above, we assume that the nonthermal electron dis- 
tribution has the horseshoe-like shape (the distributions observed 
in the terrestrial magnetosphere are shown, e.g., at Figs. 3 and 5 
in the paper of Ergun et al. 120001 1. Instead of a detailed investiga- 
tion of the formation process of a horseshoe-like distribution, we 
set the initial electron distribution to the model function similar 
to the observed ones: 



/ b (p) = Aexp 



(p-Ph) 2 



L, 

exp 



P<Pc 



(3) 



where A is the normalization factor and fi = cos a. The shape 
of the distribution function is determined by such parameters 



4 



A.A. Kuznetsov: Simulation of the Electron-Cyclotron Maser Instability 




as the typical electron momentum p^, electron dispersion in 
momentum Apb, pitch-angle boundary of the loss-cone a c (or 
p. c = cosa c ), and the loss-cone boundary width Ap. c . An exam- 
ple of the distribution function (0 can be seen at Fig. [2^. For 
a c = 0, we obtain an isotropic ring-like distribution. 

The initial energy density of plasma waves is assumed to 
equal the level of thermal oscillations: 

W^(k,t)\ „ = ^T, (4) 

where ks is the Boltzmann constant and Tq is the effective 
plasma temperature. In addition, it is assumed that the energy 
density of plasma waves cannot fall below the thermal level 
(HJl during the process of wave/particle evolution, due to spon- 
taneous radiation. 

In all simulations, we assume that the nonthermal distri- 
bution function (0 has Ap^/p^ - 0.2 and Ap c = 0.2. The 
thermal component of the plasma (if present) is described by 
a maxwellian distribution with temperature of 10 6 K. The ini- 
tial temperature of plasma waves To equals 10 6 K. The remain- 
ing parameters of the considered simulation models are given 
in Table [1] they were chosen in order to explore the influence 
of various factors on the process of wave/particle evolution. In 
all cases, the plasma density is relatively low, so that the ratio 
u) v /u>b varies in the range from 10~ 3 to 10 _1 ; the total electron 
concentration n is calculated using the plasma frequency co p . The 
relative concentration of the energetic electrons n\,/n varies from 
10~ 4 tol. 

4. Results 

4.1. Relaxation of the horseshoe-like electron distribution 

4.1 .1 . The case when a thermal plasma component is absent 

Firstly, we consider in detail an example of coevolution of 
the electron distribution and plasma waves. We assume here 
that a low-temperature thermal component is absent (rib = «)■ 
Investigations of the dispersion relations for weakly-relativistic 
plasmas (Robinson 19861 |198Tb have shown that the wave dis- 
persion (for the waves propagating across the magnetic field, 
| cos 6*| <k 1) becomes similar to that in vacuum if the typical 
electron speed v e satisfies the condition (v e /c) > (u> p /u>b) 2 . Such 
a requirement is satisfied, e.g., for the particle energy E\, > 3 
keV and <jj p /wb S= 0.1. Thus we assume that the wave refrac- 
tion index equals unity both for the ordinary and extraordinary 
modes. Also we assume that the waves are elliptically polarized 



with the axial ratio of the polarization ellipse Te = cos 9 for the 
extraordinary mode and To = - 1 / cos 8 for the ordinary mode 
(TeTo = -1). The above relations follow from the magnetoionic 
theory when wg/w -> 1 and cj p /cl) — » (Melrose & Dulk 1991 ; 
Willes, Melrose, & Robinson 1994); however, we found that the 
simulation results are not very sensitive to the exact value of the 
axial ratio T a provided that for the quasi-transversal propagation 
\T E \ <K 1 (and, accordingly, |7b| » 1). 

As an illustration, the following parameters were chosen: 
magnetic field B = 1430 G that corresponds to the electron cy- 
clotron frequency of /b =4 GHz (a typical value for the radio 
emission of ultracool dwarfs), plasma to cyclotron frequency ra- 
tio Wp/wB = 10~ 3 that corresponds to the electron concentration 
n — «b = 2 x 10 5 crrT 3 , typical energy of the energetic electrons 
£b = 10 keV, and the loss-cone boundary a c = 60°. In Table [T] 
this parameter set corresponds to the model 15. 

Figure Q] shows the contours of growth rates of the ordinary 
and extraordinary modes at the initial moment (t = 0). Only the 
positive growth rates are shown; the contour levels are evenly 
distributed between zero and the maximal growth rate value for 
a given mode which is also shown at the figure. One can see 
that the most effective wave amplification takes place slightly be- 
low the fundamental cyclotron frequency. Both emission modes 
are generated mainly in the perpendicular direction with respect 
to the magnetic field, with a slight asymmetry due to the loss- 
cone feature. Growth rates decrease rapidly with an increasing 
frequency, but the wave amplification (in an oblique direction) 
can occur even above the cyclotron frequency; in this region, 
the emission directivity patterns become essentially asymmetric. 
One can also note that the maximal growth rate of the extraor- 
dinary mode exceeds that of the ordinary mode by more than 
two orders of magnitude. Both for the ordinary and extraordi- 
nary modes, there are also amplification regions near higher har- 
monics of the cyclotron frequency, but with much lower growth 
rates (they are not shown at the figure). Thus the fundamental 
extraordinary mode is strongly dominating. 

Figure [2] shows the electron distribution function at differ- 
ent times (the chosen time moments correspond to the differ- 
ent stages of the relaxation process, see comment to Fig. |4j). For 
some time after the beginning of the simulation, while the wave 
energy density is low, the distribution function remains very sim- 
ilar to the initial one (Fig.|2^). When the waves are amplified to a 
certain critical level, an effective diffusion of electrons on these 
waves begins. As a result, the electrons drift along the trajecto- 
ries of p z = const towards lower values of p ± , thus losing energy. 
At the early relaxation stage (Fig.[2j3), this results in formation of 
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-0.2 -0.1 0.0 0.1 0.2 

Pz/{ m e c ) 



Fig. 2. Time evolution of the electron distribution for the model 
15 (TableQ). 



a plateau slightly below the maximum of the initial distribution 
function (with respect to the momentum) and around a - 90° 
(with respect to the pitch angle). On the other hand, at the upper 
and lower boundaries of this plateau (with respect to the momen- 
tum), the gradient of the distribution function increases. At the 
middle relaxation stage (Fig. [2};), the "plateau" expands with re- 
spect both to the momentum and pitch angle, and we can see that 
this feature is not really flat: while the derivative df/dp± in this 
area approaches zero, the derivative df jdp z remains nonzero, so 
the distribution function gradually decreases with increasing p z . 
Nevertheless, since the maser amplification or damping of waves 
depends mainly on the derivative with respect to the transversal 
component of the momentum, a contribution of this region of 
the distribution function into growth rate is close to zero. At the 
upper boundary (with respect to the momentum), the "plateau" 
extends beyond the circle p - Pb (where the initial distribution 
function had a maximum) and the positive slope in p ± disap- 
pears; at the lower boundary, this slope is still quite large. At 
the late relaxation stage (Fig. |2jl), the "plateau" expands even 
further, so now we have df/dp^ < almost everywhere, and 
further amplification of waves nearly ceases. A positive slope 
remains only in the low-energy region, but the resulting growth 




0.975 0.980 0.985 0.990 0.995 
a/a B 



Fig. 3. Time evolution of the extraordinary waves for the model 
15 (Table The contour levels are 0.0075, 0.015, 0.03, 0.06, 
0. 125, 0.25, 0.5, and 0.9 of the maximum wave energy density. 



rate of the waves is very low. Further simulation has not revealed 
any qualitatively new features: the "hole" around p — slowly 
shrinks and the distribution function approaches asymptotically 
a stationary (saturated) state, while its rate of change decreases 
with time. 

Figure [3] shows the distribution of the extraordinary mode 
energy density in (a>, #)-space at different times; in addition, the 
maximal values of the amplification coefficient A are shown (this 
coefficient equals the ratio of the energy density at a given time 
to its initial level). The letters identifying the panels correspond 
to those in Fig. [2] but panel (a) is omitted since the wave en- 
ergy at the initial moment is negligible. For some time after the 
beginning of the simulation, the growth rate remains nearly con- 
stant, so the wave energy grows exponentially with time. As a 
result, at the onset of relaxation of the electron beam (Fig. |3j?), 
the waves are concentrated in a relatively narrow range with re- 
spect both to the frequency and the propagation direction; this 
region is much more narrow than the initial region of positive 
growth rate. When the logarithm of the amplification coefficient 
In A reaches values of about 24-25, the waves begin to modify 
the electron distribution. This results in a sharp decrease (down 
to zero) of the growth rate in those regions where it was large 
at the beginning of the simulation, and also in an increase of 
the growth rate in the adjacent regions. These changes reflect 
formation of a plateau on the electron distribution function (see 
Fig. |2j3). As a result, further increase of the maximal intensity 
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Fig. 4. Time profiles of the total wave energy (a) and average 
growth rate (b). At panel (a), solid line is the result of numer- 
ical simulations and dashed line is a simplified functional fit. 
Simulation parameters correspond to the model 15 (Table[TJ. 



the time profile of the total wave energy can be described with a 
good accuracy by the well-known saturation curve: 

W(t) = W ss + (Woo - W ss ) [l - e -^ /T -'] , (5) 

where W ss = W(t ss ). The distributions of the electrons and waves 
at Figs.[2]|3]correspond to the times t — 0, t — f ss (early stage of 
relaxation), t = t ss + T sat (middle stage), and t = t ss + 3r sat (late 
stage). For the simulation parameters used in this section, re- 
laxation begins at f ons ^ 6.5 /is or 22.6y max , where y max is the 
maximal growth rate of the extraordinary mode at the initial mo- 
ment; the time of the steepest slope f ss ^ 9.45 fis or 32.3y max ; and 
the relaxation timescale r sat 14.4 fis or 50.2y max . The final (at 
t — > oo) total energy density of the extraordinary waves equals 
Woo 4.64 x 10" 4 erg cm -3 , which amounts to 13.3% of the 
energy density of the accelerated particles at the initial moment. 
The obtained timescales are similar (by order of magnitude) to 
those given in the article of Aschwanden ( 1990). However, the 
conversion efficiency of the particle energy into waves is now 
considerably higher, which is caused by the different type of un- 
stable electron distribution (Aschwanden ( 1990) considered only 
the loss-cone). 

As said before, growth rate of the ordinary mode is much 
lower than that of the extraordinary mode. Simulations consid- 
ering both emission modes simultaneously have shown that the 
ordinary mode is amplified by less than a factor of 1.25 in com- 
parison with the level of thermal oscillations. Thus the resulting 
energy of the ordinary waves is negligible and they neither make 
a measurable contribution to the radio emission of planets and 
UCDs nor affect the electron distribution. 



of waves nearly ceases, but the wave distribution in phase space 
broadens (see Fig. [3J; which corresponds to the middle relax- 
ation stage). At a later stage (Fig.[3jl), the distribution of plasma 
waves becomes even more broad. Note that the shape of the fi- 
nal wave distribution differs considerably from the shape of the 
initial region of positive growth; in particular, relaxation of the 
electron beam allows waves to be generated at noticeably lower 
frequencies than at the initial moment. On the other hand, the 
wave energy density in the frequency range above the cyclotron 
frequency remains negligible throughout the relaxation process 
although the growth rate can be initially positive in this region. 

Figure|4]shows the time history of the integral characteristics 
of the extraordinary waves: total energy density (integrated over 
the frequency and propagation angle) and the averaged growth 
rate (only the positive values of the growth rate were consid- 
ered). In general, the time history is very similar to that obtained 
by Aschwanden (1990). For some time after the beginning of 
the simulation (0 < t < f ons ), the growth rate is constant and 
the wave energy grows exponentially but still remains very low. 
Then, when the wave energy reaches a certain critical level (at 
t — fons), relaxation of the unstable electron distribution begins 
and the growth rate starts to decrease. Note that this critical level 
is considerably lower than the final wave energy. After the on- 
set of relaxation (at f ons < t < f ss ), the total energy of waves 
continues to grow with an increasing rate, but the correspond- 
ing curve somewhat differs from an exponential one. At time f ss 
(which is defined as the time of the steepest slope), the rate of 
change of the total wave energy dW/dt reaches its maximum. 
Later (at t > f ss ), both the average growth rate and the rate of 
change of the total wave energy gradually decrease, approach- 
ing zero asymptotically; at the same time, the total wave energy 
goes asymptotically to the saturation level (Woo)- At this stage, 



4.1 .2. The case when a thermal plasma component is 
present 

Now we investigate the case when a thermal plasma component 
dominates. We use the same parameters of the magnetic field 
and energetic particles as in the previous section (/b = 4 GHz, 
n b = 2 x 10 5 cirT 3 , E b = 10 keV, a c = 60°), but now a thermal 
plasma with the concentration of «o = 2 x 10 7 cm~ 3 is present 
(«b/« = 10~ 2 ), so that the plasma to cyclotron frequency ratio 
is Wp/wB = 10~ 2 . In Table [T] this parameter set corresponds to 
the model 9. For the wave parameters, we use the cold plasma 
dispersion relation (see AppendixlAl. 

Figure [5] shows the initial growth rates (at t = 0). In a cold 
magnetized plasma, the extraordinary mode is split into two 
branches: the fast (X) mode exists above the cutoff frequency u> c , 
while the slow (Z) mode exists below the resonance frequency 
aj r ((jj c > a> r ). Under the conditions considered, both the cut- 
off and resonance frequencies almost coincide with the electron 
cyclotron frequency. However, the frequency gap between the 
amplification regions of the X- and Z-modes is much larger than 
<d c - a> r since near the cyclotron frequency the waves are heavily 
damped by the thermal plasma component. The Z-mode is gen- 
erated across the magnetic field, while the X-mode is generated 
in an oblique direction (at 75°). The maximal growth rate of 
the Z-mode exceeds that of the X-mode by more than an order 
of magnitude. The ordinary (O) mode has the smallest growth 
rate; its dispersion curve is continuous, but the amplification re- 
gion is split into two by the absorption band near the cyclotron 
frequency. 

Simulations have shown that the time evolution of the elec- 
tron distribution is very similar to that for the case without a 
thermal plasma (see Fig.|2]i. Similarly, the time evolution of the 
Z-mode (which is the dominant mode) does not differ much 



A.A. Kuznetsov: Simulation of the Electron-Cyclotron Maser Instability 



7 




Fig. 5. Initial growth rates of the extraordinary (a) and ordinary (b) waves in the presence of a thermal electron component. 
Simulation parameters correspond to the model 9 (Table [TJ, and the cold plasma dispersion relation is used. 
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density (b) at the late relaxation stage. The thermal component 
of the distribution function is shown only partially. Simulation 
parameters correspond to the model 9 (TableQ]). 



from the evolution of the extraordinary mode shown in Fig. [3] 
Therefore we do not present here a complete time history of the 
waves and particles. Figure [6] shows the late relaxation stage. 
One can see that the "plateau" on the electron distribution func- 
tion (formed by the nonthermal electrons) merges with the ther- 
mal electron population at this stage. The conversion efficiency 
of the particle energy into the Z-mode waves (only the nonther- 
mal component is considered) is about 11.8%; this is slightly 
less that for the model without the thermal plasma. The relax- 
ation timescales are very similar to the values found in the previ- 
ous section. We can conclude that in a low-density plasma (with 
Wj/wb <s 1), the exact form of the dispersion relations is not 
very important, since both the vacuum-like and cold plasma dis- 
persion relations provide almost the same results. 

We can see that the generation of the Z-mode waves is very 
effective. However, these waves cannot escape from the genera- 
tion region. Thus the question is of interest: can the intensities 



of the freely propagating modes (X- and/or O-mode) reach suf- 
ficiently high levels despite of a lower growth rate? We simu- 
lated the simultaneous evolution of the different modes. Figure|7] 
shows the distribution of the X-mode energy density at different 
times. The time history is different from that of the Z-mode since 
the X-mode generation is caused by the loss-cone feature. Since 
the loss-cone is destroyed during relaxation (due to diffusion of 
particles on the Z-mode waves), the X-mode amplification in 
some regions of phase space is replaced by absorption, so that 
the region occupied by the waves shrinks with time. Simulations 
for times later than those displayed in the figure have shown 
that the maximal intensity of waves also can somewhat decrease. 
However, the main result for the X-mode is that its intensity is 
much lower than that of the Z-mode: the parameter In A does 
not exceed four, so the wave energy density exceeds the thermal 
level by not more than a factor of 50. Therefore the X-mode has 
no noticeable effect on relaxation of the electron distribution. 

Figure [8] shows the average growth rates of the Z- and X- 
modes. One can see that the growth rate of the X-mode starts to 
decrease simultaneously with that of the Z-mode. In general, the 
ratio of growth rates remains nearly the same throughout the re- 
laxation process (despite a slight secondary increase of the aver- 
age growth rate of the X-mode at the middle relaxation stage), so 
that the growth rate of the X-mode remains considerably lower 
than that of the Z-mode. We have found that the total energy den- 
sity of the X-mode waves does not exceed 10~ 13 erg cm -3 , which 
amounts to 3 x 10 _1 1 of the energy density of the accelerated par- 
ticles. The energy density of the O-mode is much lower than that 
of the X-mode. Thus, in this simulation, relaxation of the unsta- 
ble electron distribution is caused entirely by the Z-mode, and 
the intensities of the other wave modes are negligible. 



4.2. Comparison of the results for the different simulation 
models 

Simulations of coevolution of the electron distributions and 
plasma waves were performed for the various parameter sets (see 
TableQ}. The cold plasma and vacuum dispersion relations were 
used for the models with n\,/n <k 1 and n^jn = 1, respectively. 
In all cases, the extraordinary mode (or its slow branch, for the 
models including a thermal plasma) dominated strongly. TableQ] 
also contains the maximal initial growth rates, the timescales of 
the relaxation process, and the transformation coefficient of the 
particle energy into the waves. Figures l9lfT3l illustrate the effect 
of the various factors. 
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Fig. 10. Same as in Fig. |9j for the different values of the total 
plasma density (concentration of the accelerated electrons rib 
is constant). The values correspond to the models 15, 9, and 6 
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Fig. 8. Time profiles of the average growth rates of the Z- and 
X-modes for the model 9 (Table Q]). For the X-mode, the scale 
is different from that for the Z-mode to make the growth rate 
variations more visible. 



£b and the loss-cone boundary a c equal 10 keV and 60°, re- 
spectively. One can see that the relaxation timescales are simply 
inversely proportional to the concentration of the energetic par- 
ticles. The conversion efficiency of the particle energy into the 
waves increases slightly with increasing n^/n and varies from 
11.2% to 13.6%. 

Figure[10]shows the simulation results for the case when the 
concentration of the energetic particles is constant («b = 2 x 10 5 
cm -3 ), while the total plasma density varies (which results in 
different ratios both of n^/n and <jj p /wb)- With increasing w p /wb 
from 10 3 to 10" 2 , the growth rate of the extraordinary mode and 
the relaxation timescales remain nearly unchanged. At the same 
time, with increasing w p /wb from 10" 2 to 10" 1 , the growth rate 
decreases by about an order of magnitude (due to the chang- 
ing dispersion parameters of the waves), while the relaxation 
timescales increase by the same factor. The conversion efficiency 
of particle energy into waves varies around 12% and slightly de- 
creases with increasing w p /wb- 



4.2.1 . Effect of varying the plasma parameters 

Figure [9] shows the simulation results (relaxation timescales and 
conversion efficiency of the particle energy into the waves) for 
the different relative concentrations of the energetic electrons 
(rib/n = 10" 4 , 10" 2 , and 1) while the total plasma density is 
assumed to be constant (n = 2 x 10 7 cm" 3 that corresponds to 
u) p /u>b = 10" 2 ). The typical energy of the accelerated electrons 



4.2.2. Effect of varying the distribution of the energetic 
electrons 

In this section, we consider different distributions of the ener- 
getic electrons. In all cases, the concentration of the energetic 
electrons is assumed to be the same («b = 2 x 10 5 cm" 3 ). Both 
the models without a thermal plasma component and the mod- 
els including such a component (with «b/n = 10" 2 ) are consid- 
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Fig. 12. Same as in Fig. [9] for the different angular distributions 
of the accelerated electrons. Solid lines: models 14, 15, 16, and 
17 (without thermal component); dashed lines: models 8, 9, 10, 
and 1 1 (with thermal component). 



ered; the corresponding points at Figs. [TTIfT2l are connected by 
the solid and dashed lines, respectively. 

Figure QT| shows the simulation results for the different typ- 
ical energies of the accelerated particles (in all cases, the loss- 
cone boundary is a c = 60°). For the models without a thermal 
plasma component, the conversion efficiency of the particle en- 
ergy into the waves is almost constant (about 13.3%). For the 
models including a thermal plasma, the conversion efficiency in- 
creases with E b (from 8% at E h = 3 keV to 13% at E h = 30 
keV). An increase of the electron energy makes the relaxation 
process slower (both the timescales f ss and r sat increase). This 
is because the "plateau" formation for the electron distributions 
with higher energy requires a larger displacement of particles in 
the momentum space and, consequently, a higher energy density 
of plasma waves (in order to provide a stronger diffusion) and a 
longer time. For the models with a thermal component, the re- 
laxation process is slower. With the increasing electron energy, 
the relaxation timescales (as well as the energy conversion effi- 
ciency) for the models with a thermal component approach the 
corresponding values for the models without a thermal compo- 
nent. 

Figure[T2]shows the simulation results for the different loss- 
cone boundaries (in all cases the beam energy is E b = 10 keV). 
The models without a thermal plasma component always pro- 
vide a higher conversion efficiency of the particle energy into the 
waves than the corresponding models with a thermal component. 
The highest conversion efficiency as well as the fastest relaxation 
occur for a c = 60°. For the ring-like distribution (with or c = 0), 



as well as for the loss-cone with a c = 90°, the conversion ef- 
ficiency and relaxation rate are slightly lower. The distribution 
with a c = 120° is similar to those used in the paper of Bingham 
& Cairns d2000l > and Bingham, Cairns, & Kellett d200TT) : it has 
no electrons with pitch angles around 90°. Therefore, for this 
distribution, the amount of free energy is relatively low and only 
about 5-6% of the initial electron energy can be transferred to 
waves. We would like to highlight that for all considered distri- 
butions (including essentially anisotropic ones), the waves are 
generated preferably in the perpendicular direction to the mag- 
netic field; e.g., for the distribution with a c = 120°, the max- 
imum of the wave intensity (at the late relaxation stage) is at 
6*^91°. 

4.2.3. Effect of varying the magnetic field 

In this section, we investigate the influence of the magnetic 
field strength on the coevolution of the electron distributions and 
plasma waves. Figure [T"3l shows the simulation results for three 
values of the electron cyclotron frequency: 400 kHz (which is 
typical for the AKR sources in the terrestrial magnetosphere), 
40 MHz (which is typical for the magnetosphere of Jupiter), and 
4 GHz (which seems to be typical for the magnetospheres of 
UCDs). We assume that the plasma to cyclotron frequency ratio 
is constant and equals u> p /cl>b = 10~ 3 for the models without a 
thermal plasma component; thus the concentrations of the ener- 
getic electrons rt\, equal 2xl0~ 3 , 2X10 1 , and 2xl0 5 crrT 3 , respec- 
tively. In the models including a thermal plasma, we assume that 
« p /«b = 10 -2 and «b/« = 10~ 2 which provides the same con- 
centrations of the energetic electrons as above. The typical en- 
ergy of the accelerated electrons E\, and the loss-cone boundary 
a c equal 10 keV and 60°, respectively. Under the conditions con- 
sidered, the growth rates of the plasma waves are proportional to 
the cyclotron frequency (see TableQ]). On the other hand, the re- 
laxation timescales decrease with an increase of the cyclotron 
frequency somewhat faster than the inverse proportionality law 
implies. The largest relaxation timescales (relative to the inverse 
cyclotron frequency) occur at the Earth since in a weaker mag- 
netic field diffusion of particles on plasma waves also weakens, 
and the relaxation onset requires a higher wave amplification co- 
efficient A (this effect is equivalent to that of reducing the initial 
temperature of plasma waves). The fastest relaxation (both in 
absolute and relative units) occurs for UCDs. A complete re- 
laxation of the unstable electron distribution is achieved in less 
than 1 s (at the Earth), less than 10 ms (at Jupiter), and less than 
0.1 ms (for UCDs). The conversion efficiency of the particle en- 
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ergy into the waves slightly increases with an increasing mag- 
netic field strength; the models without a thermal plasma compo- 
nent always provide a higher conversion efficiency (13.1-13.3%) 
than the corresponding models with a thermal component (1 1.8- 
12.0%). 

5. Discussion 

Our simulations have shown that in a relatively low-density 
plasma with a> p /a)s = 10~ 3 - 10 _I , the ring-like or horseshoe- 
like distributions of accelerated electrons generate mainly the 
extraordinary waves with the frequency slightly below the elec- 
tron cyclotron frequency. This conclusion is consistent with 
the results of previous studies (e.g., Ergun et al. I20001 I. If the 
magnetoionic theory is applicable (i.e., the wave dispersion is 
dominated by a cold plasma component), the generated waves 
correspond to the Z-mode branch of the extraordinary mode. 
Generation of other modes (such as the ordinary mode, the X- 
mode of the magnetoionic theory, and the waves near higher cy- 
clotron harmonics) is also possible, but their final energy den- 
sity is less than that of the fundamental extraordinary mode by 
many orders of magnitude, even if the difference in the initial 
growth rates is not so large. We have found that the dominat- 
ing mode remains the same throughout the relaxation process of 
an unstable electron distribution. This differs from the results of 
Fleishman & Arzner (2000), where relaxation of a loss-cone dis- 
tribution on the initially dominating mode (lower-hybrid waves) 
resulted in formation of a distribution which was stable with re- 
spect to excitation of lower-hybrid waves but still sufficiently 
anisotropic to amplify other types of waves (X and O). Most 
probably, this is because Fleishman & Arzner (2000) consid- 
ered a different electron distribution (loss-cone vs. horseshoe) 
and different plasma parameters (u) p /(l>b > 0.2). In our simu- 
lations, relaxation of an electron distribution on the waves of 
the dominating mode prevented amplification of other modes as 
well. Also note that for other types of unstable electron distribu- 
tions, another wave mode can prevail; for example, a loss-cone 
distribution in a low-density plasma excites mainly the oblique 
X-mode (Aschwanden |1990j l. 

For the ring-like or horseshoe-like electron distributions, the 
transformation coefficient of the particle energy into waves can 
exceed 10%. This is considerably higher than for the loss-cone 
distributions (Aschwanden 1990). The reason is that for the loss- 
cone, only a small fraction of particles (with pitch-angles around 
a c ) actually makes a contribution to the generation of waves and 
provides them with energy; in contrast, for the horseshoe-like 
distribution, the process of wave generation involves the major- 
ity of particles. It is interesting to note that the conversion effi- 
ciency is very weakly dependent on the plasma parameters and 
is determined mainly by the shape of the nonthermal distribu- 
tion. The presence of a thermal plasma component reduces the 
conversion efficiency. 

We estimate now the required efficiency of the emission 
mechanism. As an example, we consider the M9 dwarf TVLM 
513^-6546 located at a distance of d = 10.5 pc. This object is 
known to produce radio bursts with an intensity up to / — 6 mJy 
(at the frequency about 8 GHz) (Hallinan et al. 2007). The typ- 
ical size of the AKR generation region (in the direction across 
the magnetic field) can be estimated as 300 km. The radius of 
the UCD is expected to be comparable with that of Jupiter (i.e., 
ten times larger than the Earth), so we can assume that the radia- 
tion source has a size about R ± 3000 km. Thus the brightness 
temperature of the emission is about 10 13 K. If the initial temper- 
ature of plasma waves equals To = 10 6 K then the logarithm of 



the amplification coefficient should be not less than In A 16. 
For the growth rate y 3 x 10 6 s _1 (see Section 14.1.1 1 for the 
beam and plasma parameters) and group velocity of the waves 
u gr < c, such amplification can be achieved at a distance of no 
more than 1 .6 km. Thus the estimated growth rates of the ECMI 
are well in excess of the required values, even for the relatively 
low concentrations of the accelerated electrons. 

On the other hand, as stated above, estimations based only 
on the growth rate are insufficient and we need to check the con- 
version efficiency of the electron energy into radiation. The total 
power of radio emission from the UCD (for the case of isotropic 
radiation) can be estimated as F r = 4nd 2 IAf, where A/ is the 
spectral bandwidth of the emission. For the above parameters, 
assuming that Af / 8 GHz (this is an upper limit as the 
spectrum is unlikely to be this broad), we obtain F, - 6x 10 21 erg 
s . The energy flux of the accelerated electrons can be estimated 
as fb = n\,v\,E\,R\, where «b, v\>, and E\, are the concentration, 
speed, and energy of particles, and R\ is the cross-section area 
of the electron beam. We assume that the plasma to cyclotron 
frequency ratio is u> p /u>b = 10~ 3 (which is less than that in the 
AKR sources) and all particles in the radiation source are accel- 
erated («b = n). Then for the cyclotron frequencies f% =4-8 
GHz we obtain «b — (2 - 8) x 10 5 cm 3 . For the electron en- 
ergy £b = 10 keV, the corresponding speed i>b = 5 x 10 9 cm 
s _1 , and the beam size R ± = 3000 km, the energy flux of the 
particles is F\, (1.4 - 5.7) x 10 24 erg s . Thus the conversion 
efficiency of the electron energy into the radio emission should 
be not less than F r /Fb - (1 - 4) X 10~ 3 . According to our sim- 
ulations, the ECMI efficiency (> 0.1) is much larger. The above 
estimations do not take into account the possible absorption of 
the radio emission during propagation, radiation directivity, and 
(possibly) a narrow spectral band. Nevertheless, with reasonable 
assumptions about the source parameters, the ECMI is well able 
to provide the observed intensity of radio emission from UCDs. 

The terrestrial AKR is generated in the auroral cavity where 
a cold plasma component is almost absent; this allows the waves 
produced below the electron cyclotron frequency to escape di- 
rectly from the source region due to relativistic corrections to 
the dispersion relation. It is very likely that the radio emis- 
sion of UCDs is generated under similar conditions. However, 
let us discuss a possible case when a cold plasma component 
dominates. Under such conditions, the waves (of the Z-mode) 
excited by the ECMI cannot escape from the source due to a 
stop band at a> - a>s- Thus the Z-mode has to be transformed 
into electromagnetic radiation, e.g., due to nonlinear processes. 
We can estimate the required efficiency of nonlinear conversion 
as 77nl ^ (1 - 4) x 10" 2 (for the case when w p /«b = 10~ 2 , 
«b/« = 10~ 2 , 10% of the particle energy goes into the Z-mode 
waves, and all other parameters are the same as in the previous 
paragraph). Nonlinear processes require a further investigation, 
but at present we cannot rule out the described two-stage emis- 
sion mechanism. It is interesting to note that the growth rates of 
the X- and O-modes (see Fig. [5), in a linear approximation, also 
can provide the required amplification coefficient at distances of 
tens to hundreds of kilometers. However, the kinetic simulations 
have shown that the conversion coefficient of the particle energy 
into these waves is extremely small, so a direct generation of 
electromagnetic waves by the horseshoe-like electron distribu- 
tion in the presence of a cold plasma cannot be responsible for 
the observed emission. 

According to the simulation results (see Table [TJ, the typi- 
cal relaxation time of electron distributions with w p /wb = 10~ 2 
and rib In a 1 in the terrestrial magnetosphere should be about 
^diff — T"sat — 2 x 10~ 3 s. On the other hand, the typical escape 
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time of the waves from the generation region is r esc /?j_/t> gr . 
For R ± a 300 km and u gr a c, we obtain r esc - 10~ 3 s. Thus 
Tdift > fesc, and the diffusive limit is, in fact, inapplicable. In 
the terrestrial magnetosphere, escape of waves from the ampli- 
fication region has a dominating effect on the relaxation of the 
electron beams. In particular, this allows the strongly unstable 
electron distributions with large values of df/dp± (such as re- 
ported by Ergun et al. [2000) to exist in a quasi-stationary state. 
Escape of waves from the amplification region obviously re- 
duces the relaxation rate, so that the processes of particle acceler- 
ation and magnetic mirroring are able to compensate the changes 
in the distribution function caused by diffusion on plasma waves. 
On the other hand, for UCDs we obtain r esc a 10~ 2 s, while 
Tdiff - 10~ 7 - 10~ 5 s; thus Tdiff <k T esc , and the use of a dif- 
fusive limit for UCDs is well justified. We can expect that the 
quasi-stationary electron distributions in the magnetospheres of 
UCDs will be very close to the relaxed state (as shown, e.g., 
at Figs. |2j;-d). To get more accurate results, one has to include 
the mechanism responsible for formation of an unstable electron 
distribution into the simulation model. The case when escape of 
waves from the amplification region is important for the relax- 
ation process will be considered in a future paper. 

Relaxation of unstable electron distributions in the magne- 
tospheres of UCDs is very fast (<SC 1 s). On the other hand, the 
observed bursts of radio emission have a much longer duration 
(tens of seconds). It is very likely that the emission source is ac- 
tually even more long-lived, but we observe the bursts only when 
the narrowly-directed radio beam (whose direction changes due 
to rotation of a star) points towards the observer (Hallinan et 
al. 120061 120071 Berger et al. [2009l l. This again means that the 
diffusive limit can be considered only as a very rough approxi- 
mation. Actually, the emission properties are determined by the 
relatively stationary (but long-living) processes of particle accel- 
eration; any model attempting to make a quantitative interpreta- 
tion of observations has to include these processes. Nevertheless, 
we believe that the results obtained in this paper allow us to 
make some conclusions about the processes in magnetospheres 
of UCDs, and will be useful in developing more advanced mod- 
els. 

6. Conclusion 

In this paper, we made a numerical simulation of the electron- 
cyclotron maser instability and investigated coevolution of un- 
stable electron distribution and plasma waves in a low-density 
plasma (w p /wb <k 1). Electron distribution of the "horseshoe" 
type was considered as an energy source for the plasma oscilla- 
tions. The simulation parameters were chosen according to the 
measurements in the terrestrial magnetosphere, and scaled to the 
conditions expected in the magnetospheres of UCDs (where the 
electron cyclotron frequency is a few GHz). Spatial movement 
of the waves and particles, as well as other processes modifying 
the electron distribution (except of diffusion on plasma waves), 
were neglected. A number of simulations were made for the dif- 
ferent parameters of plasma and accelerated particles. We found 
that: 

- Under the conditions considered, the electron beams with 
ring-like or horseshoe-like distributions generate mainly the 
extraordinary waves propagating across the magnetic field, 
with a frequency slightly below the electron cyclotron fre- 
quency. Other wave modes also can be amplified but their 
final intensity is less than that of the fundamental extraordi- 
nary mode by several orders of magnitude. 
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- Conversion efficiency of the energy of accelerated electrons 
into the extraordinary waves can exceed 10%; this parameter 
depends mainly on the shape of the nonthermal distribution. 

- A typical relaxation time of the unstable electron distribution 
is a few tens of inverse growth rates. In the magnetospheres 
of UCDs, a complete relaxation of the electron distribution 
should be very fast - much less than one second; thus the ob- 
served emission requires a long-lived source of accelerated 
particles. 

- Energy estimations show that the efficiency of the electron- 
cyclotron maser instability is sufficient to provide the ob- 
served intensity of radio emission from UCDs under reason- 
able assumptions about the energies and concentrations of 
accelerated particles. 
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Appendix A: Dispersion of the magnetoionic modes 

The refraction index of the electromagnetic waves (No-) in a cold 
magnetized plasma satisfies the dispersion equation 



kc 

Ni = \ — 
where 



1 



2V(1 - V) 



2(1 -V)- U sin 2 O + o- V© 



(A.l) 



By expanding the polarization parameters dA.411A.5l l in the small 
parameters yfU - 1 and y/V and retaining only the zero or- 
der terms in the expansions, we obtain the polarization state of 
the magnetioinic modes near the cyclotron frequency in a low- 
density plasma: 

r_i - cos0, T + i =* — , L ±x * 0. (A. 12) 

cos 



D=U 2 sin 4 + 4U(1 - V) 2 cos 2 6, 



U = (« B /«) 2 , V = (co p /to) 2 , 



(A.2) 
(A.3) 



to and k are the wave frequency and wave vector, a) p and wb are 
the electron plasma and cyclotron frequencies respectively, and 
is the angle between the wave vector and the magnetic field. 
In Eq. (1A.1I ). cr = -1 for the X-mode (fast extraordinary) and 
the Z-mode (slow extraordinary); cr = +1 for the O-mode (fast 
ordinary) and the W-mode (whistlers). 

The polarization state of the magnetoionic modes is de- 
scribed by the following parameters: 



Lrr = 



2 Vt7(l - V)cos6» 
U sin 2 9 - cr V© ' 

V Vt7sin0+ 7Vt/Vsin0cos0 
1 - U -V + UVcos 2 ' 



(A.4) 



(A.5) 



where T a is the axial ratio of the polarization ellipse and L a is 
the longitudinal part of the polarization. 

The group velocity of the magnetoionic waves equals 



..("■) 



da) c 
" gr ~dk ~ d{a)N lT )lda) 

where 



dto 



1 + 



V yfUT a cos 



2(7V- Vf/costf) 2 

(i + y)d-r 2 )" 



(1 - V)(l + Tl) 



(A.l) 



The derivative of the refraction index with respect to the propa- 
gation direction can be calculated using the relation 



1 dN a 



1 T 

^(T 1 {J 



80 l + ri 



(A.8) 



The different magnetoionic modes exist in the following fre- 
quency ranges: 

- X-mode: at u> > u> c+ ; 

- O-mode: at u> > u> p ; 

- Z-mode: at <w c _ < a> < a> r+ ; 

- whistlers: at a> < w r _. 

The cutoff and resonance frequencies are given by: 



1 

OJ c± = ±2 W b + 



W r± = 2^ W B + w p) ± y 4(^6 + W p) ~ W B W P C0S ' 

In vacuum, the dispersion parameters become: 



No- = 1, 



„("■) 



d(ojN IT ) 
da> 



8N a 
80 



(A.9) 



(A. 10) 



(A.ll) 



Appendix B: Expressions for the growth rate 

According to Melrose & Dulk dl982.l t and Aschwanden d 19901 1, 
the growth rate of the magnetoionic oscillations equals (the 
wave-mode index cr is hereafter omitted for brevity) 

47rVc 2 



r (k) = 



jNd(o)N)/dco(l +T 2 ) 



Z 



7\cos 0- N/3 : ) + L sin0 



NJ3 ± 
soj B 8f(p) 



J S (A) + J'M) 



dp- Tv ± dp ± 
d(u>-k z v z - d 3 p, 



Pi 



(B.l) 



where v and p are the electron velocity and momentum, the in- 
dices z and _L indicate the longitudinal and transversal compo- 
nents of the vectors with respect to the magnetic field (in partic- 
ular, N z = Ncos0andN ± = Nsin0),/3 = v/c, T = (1 -yS 2 r 1/2 
is the relativistic factor, J S (A) and J' S (A) are the Bessel function 
and its derivative over the argument A, 



(A.6) A 



k±p± 
m e UB 



(B.2) 



and /(p) is the electron distribution function satisfying the nor- 
malization condition 



/(P)d j p = «, 



(B.3) 



where n is the total electron concentration. The growth rates of 
the vacuum modes can be obtained by substituting the corre- 
sponding dispersion parameters (1A.11HA.T21 . 
If we introduce the dimensionless parameters 



u = — , x= — , Y = — , 



(B.4) 



and, in addition, use the polar coordinates (u, a) for the distribu- 
tion function (where a is the angle between the electron velocity 
and the magnetic field), then the expression for the growth rate 
takes the form 

nY 2 



y(k) 



Nd(a)N)ldco(\ + T 2 ) 



J 

dm 



T(cos - Nfi z ) + L sin0 



NJ3 A 



J S (A) + J'M) 



~-+(cos a -N z (3)^^- 
ou da 



since 



X 5 



xN z u z 



r 



d 3 u. 



(B.5) 



The dimensionless distribution function /(u) should satisfy the 
normalization condition 



/(u)d j u = 1, 



(B.6) 
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and A = xN ± u±. 

Using the properties of the ^-function, we can reduce three- 
dimensional integrals in ( 1B.51 ) to one-dimensional integrals over 



du z : 



2n z Y 



2 V 2 



y(k) = 

w B xNd{u)N)ldu){\ + T 2 ) 



x 



;=-co u 

U z ni ; n 

em 



7\cos 9- Nf3 z ) + L sin6» 



du 

X r sin a 



NJ3 A 



+ (cos a - N-J3) 



J S (A) + J' S (A) 



em 



da 



du z , 



(B.7) 



u ± =u ± (u z ) 



where the value u ± at every point of the resonance curve is found 
from the resonance condition: 



u]_(u z ) = ^- + N z u z \ - u\ - 1 



(B.8) 



and the integration limits u zm ; m and u zmax (which are depen- 
dent on the harmonic number) are the boundaries of the interval 
where Eq. ( IB. 81 ) has a solution. 

For |jV z | < 1 (which is always satisfied for X- and O-modes), 
the resonance curve in the momentum space ( IB. 81 ) is an ellipse 
intersecting the axis u ± — at two points: 



u zl,2 



sNJx + Jn 2 + s 2 /x 2 - 1 



l-N 2 



(B.9) 



where u z \ and u Z 2 correspond to signs "— " and "+", respectively. 
In this case, M zm ; n = u z \ and u zmax = u Z 2- For \N Z \ > 1, Eq. ( IB .St 
describes a hyperbola with only one branch being physical (for 
which, at m — > oo, the longitudinal component of the momentum 
u, has the same sign as N z ). Therefore the integration limits in 
( IB. 71 ) will be equal (note that u z \ > u Z 2 in this case): u zm [ n = — oo 
and M ;max = u z 2 at N z < 0; M, m ; n = u z \ and « zmax = oo at N z > 
0. Infinite integration limits can be avoided since the electron 
distribution function is usually defined only in a finite range of 
momentums, e.g., at u < Mhigh- In this case, the infinite limit 
(upper or lower, depending on the sign of N z ) should be replaced 
by the coordinate of the intersection point of the resonance curve 
with the circle u = Mhigh, that is 



^zhigh — 



high 



■ s/x 



N- 



(B.10) 



where Fh; g h is the relativistic factor of electrons with the momen- 
tum M high . 

The above formulae for the growth rate involve infinite sums 
over cyclotron harmonics s. Actually, only the harmonics in a 
certain range s m \ n < s < s max make a contribution to the growth 
rate, since for the other harmonics either the resonance condition 
(for a given wave parameters) is never satisfied or the resonance 
curve lies outside the domain of the electron distribution func- 
tion. In this work, the interval of harmonic numbers is taken with 
a large excess (say, from -100 to 100), and then we check for 
each harmonic whether it makes a contribution into the growth 
rate; this method is simple and very fast, and allows us to take 
into account all harmonics having an effect on the growth rate. 



Appendix C: Expressions for the quasi-linear 
diffusion rate 

The change in the electron distribution function due to diffusion 
on the magnetoionic waves is described by the last equation of 
the system ([TJ. In polar coordinates (p, a), this equation takes 
the form (Aschwanden & Benz 1988 ) 



1 



3 



dm = 

dt p 2 sin a da 
)_d_ 

p 2 dp 



D aa (p)—^— + pD ap (p) 



da 



dp 



+ 3— \P 



Dpa(p) mi +pDpp(v) d f^ 

da 



dp 



(C.l) 



where D pp , D pa , D ap , and D aa are the components of the diffu- 
sion tensor (D pa = D ap ). We note that the total diffusion tensor 
equals the sum of the diffusion tensors on separate modes (for 
simplicity, the formulae below refer only to one mode). 

If we use the dimensionless momentum u and introduce, ac- 
cording to Aschwanden & Benz (1988), generalized diffusion 
coefficients 



D, 



D„ 



(m e c) 2 ' (m e c) 2 (m e c) 2 ' 

then Eq. ( IC.U can be written in a form 



D 2 



(m e c) 2 



(C2) 



df(u, a) 
dt 



1 



df(u, a) 

' — ^ 

du 



2Dq{u, a) + D\{u, a) cot a + 



dD[(u, a) 



da 



df(u, a) 
da 



3Di(u,a) 

D\(u, a) + u h U2{u, a) cot a 



du 

~ n , , d 2 /(«, a) ^ 2 9 

ZuD\(u, a) — — — V u — 

duda du 



Dq(u, a) 



df(u, a) 



du 



d_ 

da 



D 2 (u, a) 



3f(u, a) 



da 



(C3) 



Here the expression elements are rearranged to make the numer- 
ical calculations more convenient and accurate. 

The diffusion coefficients D r (for an individual wave mode) 
are (Aschwanden & Benz ll988l Aschwanden ll990l) 



4-n e sin a v~i 
A-(u) = T1 — J] 

C S — — OQ 

W k (k) 



cos a — NJ3 ) 



Nd{u)N)ldu){\ + T 2 ) 
~T(cos6-N,3 z ) + Lsm0 



N ± /3 ± 

5 1 cd - k z v z —J a k, 



JM) + J'M) 



(C4) 



where Wk is the energy density of the considered mode in the 
space of wave vectors. The total energy density (the energy per 
unit volume) equals 



W = j W k (k)d 3 k. 



(C5) 



The diffusion coefficients for the vacuum modes can be obtained 
by substituting the corresponding dispersion parameters (I A. 111 - 
1X121 ). 
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If we use in ( lC.4t the dimensionless frequency x and the 
propagation direction as the wave characteristics instead of the 
wave vector, then the diffusion coefficients become 



7 2 • 2 00 

9 kue sin a v~i 

D r (u) = ^ 2 r 2 V 



T(cos 6 -N0 z ) + L sine 



cos a - NzfiV 



-_N£X 
a I 



NJ3 ± 
< W(x,0)Nx 2 sin 
1 + T 2 



6\x 



J S (A) + J' S (A) 
s + xN,u 



dxd6, 



(C6) 



where "W is the relative (with respect to the level of thermal os- 
cillations) energy density of plasma waves in the space of wave 
vectors: 



W(x,0) = 



W k (x,0) 



W 



w (0) = ^T 



(2n) 



3 ' 



(C7) 



According to the initial conditions used, at the initial moment 
we have 'Wix, 0) = 1 . 

The integral over d0 in (IC.6b can be calculated analytically 
by using the properties of the 5-function (see Aschwanden & 
Benz 1988; Aschwanden 1 1 990b . As a result, the expression for 
the diffusion coefficients takes the form 



Dr(u) 

a> B 



o) B T 



k-&e sin a 
m e 2 c 5 jftj 



s 



cos a - N : J3\ r 



-JM) 

a I 



7\cos - N/3 Z ) + L sin 



< W{x,0)x sin 6» 



/ S (A) + y;c^) 



(1 + T 2 )| sin - (1 /N)(dN/d0) cos 0| 



dx, (C.8) 



e=B(x) 



where the integration limits x m \ n and x max correspond to the do- 
main of the function "W(x, 0). 

The resonance curve (dependence of on x) is found using 
the following considerations. If the particle momentum u and 
wave frequency x are known, then the resonance condition re- 
sults in 



N, = 



r - s/x 



(C9) 



If the wave dispersion is assumed to be like in vacuum, then 
N - I and we immediately obtain cos — N z (the requirement 
| cos 0\ < 1 must be satisfied). If the cold plasma dispersion rela- 
tion is used, then the solution becomes a bit more complicated. 
With a known N z - N cos 0, the dispersion equation ( IA.1I ) can be 
transformed to the biquadratic equation in the variable r\ = cos 0: 



0, 



Wrf + Srj 2 + C 
where 

J[ = U - (1 - V) 3 - UV(1 - N 2 ), 

£ = N 2 [2(l - V) 2 + UV(l - N 2 ) - 2U], 

C = N*(U + V - 1). 

This equation has a solution 



2M ' 



(CIO) 

(Cll) 
(C12) 
(C.13) 

(C14) 



In order to determine the correct sign ("+" or "-") in the above 
formula for a given mode, one has to check whether the obtained 
values of the angle satisfy the dispersion equation ( IA.U . In 
practice, it is sufficient to check the validity of the inequality 



2y(1 - y) " 2(i-y ) + f/(i-, 2 ) 



T] 2 -N 2 



> 0. 



(C.15) 



In addition, the obvious requirements S 2 - AJ[C > and rj 2 < 1 
must be satisfied. Depending on the conditions, Eq. (IC.lOb for 
a given magnetoionic mode can have up to two solutions. If a 
solution with respect to rj 2 exists then the sign of 77 coincides 
with that of N z dOt. 

The interval of harmonic numbers contributing into the dif- 
fusion coefficients (s m i n < s < s max ) can be estimated analyti- 
cally, by using the resonance condition. As a result, for X- and 
O-modes (which have the refraction index N < 1) we obtain 



1 < s < 2Tx„ 



(C.16) 



For the Z-mode and whistlers, the refraction index is not limited 
from above (in the cold plasma approximation), and we have to 
introduce such limits artificially: let < AWx- Then we obtain 



(C.17) 



In this work, we use Af m ax = 10 which exceeds the actual values 
for the waves in the region of positive growth rate. The above 
formulae, as a rule, estimate the range of harmonic numbers with 
a large excess. 



Appendix D: Numeric code 

In our simulations, the electron distribution function and energy 
density of plasma waves are represented as two-dimensional ar- 
rays: 



(D.l) 



whereO < i < N x -l,0 < j < N g -l,0 <m< M,-l,and0 < n < 
N a -l.Asa result, the growth rate and diffusion coefficients need 
to be calculated only in the nodes of the corresponding grids, 
and these values (for each considered mode) can be represented 
as two-dimensional arrays as well: 



y(x, 0) -> jij, D r (u, a) 



D' 



(r) 



(D.2) 



Since the growth rate and diffusion coefficients are linearly de- 
pendent on the electron distribution function and wave energy 
density, respectively (see Eqs. JB.lt and (IC.4t ). they can be com- 
puted using tensor multiplication (Fleishman & Arzner l2~000t : 



Yij — Rijmnfmn i 



(f) 



r ijmn vy 'J> 



(D.3) 



where the kernels R; ,„,„ and P . . are constant throughout the 

ijmn t j mn to 

simulation process. Therefore they can be computed once which 
reduces the computation time considerably. 

Since the resonance curves in the spaces of particle momen- 
tums and wave vectors, as a rule, do not pass through the grid 
points, the calculation of the growth rate and diffusion coeffi- 
cients in the discretized case requires interpolating. In this work, 
we use a bilinear interpolation for the energy density of plasma 
waves. For the electron distribution function, we actually need 
to calculate its partial derivatives; this is done using a linear 
interpolation in one variable and a Lagrange interpolation on 
three closest points in another variable (where the derivative is 
needed). 
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For the electron distribution function, we use the grid where 
the momentum nodes are evenly distributed in the open inter- 
val from mi ow to Mhigh, and the pitch-angle nodes are evenly dis- 
tributed in the open interval from to n, so that 



For this grid, the best results are achieved (i.e., the conservation 
laws for the total energy and particle number are fulfilled with 
the highest accuracy) with the following boundary conditions 
(which are numerically implemented as an extrapolation of the 
arrays f m „ and D^ n beyond the grid): 

f-l,n = fo.n, h u ,n = //V„-l,n> 



fm-1 — /m,0> fm,N„ - fm,N„-\, 



-l,n 




D (0) = 

N u ,n 


D (0) 
Af„-1,«' 


ZJ (1 > 

-l,n 


0,n' 


D m = 

N u ,n 


D m 

N u -l,n' 


m,—\ 


= -D (1 l 

m,0' 




= -D ( \ 

m,N a - 


m-l 


m,0' 


m,N„ 


= ~D ( \ 

m,N a - 



(D.5) 



For the energy density of plasma waves IV,;, boundary condi- 
tions are not needed since the equations used do not contain 
derivatives of this value with respect to the wave parameters. 




(D.4) 



